% +ROMETEST\MODEL_TEST_CONDITIONAL_CONSTRAINT Test Script to model and
% solve a simple conditional constraint problem. 
%
% This model is a modification of the portfolio optimization model used in
% model_test_price_of_robustness.
%

% Begin Rome Model
disp(sprintf('\nTesting Conditional Constraints ... '));

% Paramter setup
n  = 150;                       % Number of stocks 
mu  = 1.15+ 0.05/150*(1:n)';    % Mean return
sigma = 0.05/450*sqrt(2*n*(n+1)*(1:n)'); % Deviation
Gamma = 4;  

robust_obj = zeros(3,1);
exp_return = zeros(3,1);

%% Case 1
disp(sprintf('\nCase 1 ...'));

% begin rome
h = rome_begin;   

% Bertsimas and Sim's uncertain set 
newvar z(n) uncertain;
rome_box(z, -1, 1);    

% Returns relations with uncertain factors
r = mu + sigma.*z;

% Portfolo weights
newvar x(n) nonneg;
rome_constraint(sum(x)==1); 

[zz U1] = rome_clone_uncertainty(z);
rome_constraint(norm1(zz) <= Gamma);

[zz U2] = rome_clone_uncertainty(z);
rome_constraint(norm1(zz(1:2:150)) <= Gamma);

[zz U3] = rome_clone_uncertainty(z(2:2:150));
rome_constraint(norm1(zz) <= Gamma);

% Objective
newvar s
rome_conditional_constraint(s <= r'*x, U1);

rome_maximize(s);   % Note that r is uncertain and depends in z.     
                    % The objective should be interpreted as 
                    % max{y : y<= r(z)'x for z in uncertainty set}

% solve
h.solve;
robust_obj(1) = h.objective;
exp_return(1) = mu'*h.eval(x);

rome_end;

%% Case 2
disp(sprintf('\nCase 2 ...'));

% begin rome
h = rome_begin;   

% Bertsimas and Sim's uncertain set 
newvar z(n) uncertain;
rome_box(z, -1, 1);    

% Returns relations with uncertain factors
r = mu + sigma.*z;

% Portfolo weights
newvar x(n) nonneg;
rome_constraint(sum(x)==1); 

[zz U1] = rome_clone_uncertainty(z);
rome_constraint(norm1(zz) <= Gamma);

[zz U2] = rome_clone_uncertainty(z);
rome_constraint(norm1(zz(1:2:150)) <= Gamma);

[zz U3] = rome_clone_uncertainty(z(2:2:150));
rome_constraint(norm1(zz) <= Gamma);

% Objective
newvar s
rome_conditional_constraint(s <= r'*x, U2, U3);

rome_maximize(s);   % Note that r is uncertain and depends in z.     
                    % The objective should be interpreted as 
                    % max{y : y<= r(z)'x for z in uncertainty set}

% solve
h.solve;
robust_obj(2) = h.objective;
exp_return(2) = mu'*h.eval(x);

rome_end;

%% Case 3
disp(sprintf('\nCase 3 ...'));

% begin rome
h = rome_begin;   

% Bertsimas and Sim's uncertain set 
newvar z(n) uncertain;
rome_box(z, -1, 1);    

% Returns relations with uncertain factors
r = mu + sigma.*z;

% Portfolo weights
newvar x(50) nonneg;
rome_constraint(sum(x(1:50))==1); 

[zz U1] = rome_clone_uncertainty(z);
rome_constraint(norm1(zz) <= Gamma);

[zz U2] = rome_clone_uncertainty(z);
rome_constraint(norm1(zz(1:2:150)) <= Gamma);

[zz U3] = rome_clone_uncertainty(z(2:2:150));
rome_constraint(norm1(zz) <= Gamma);

% Objective
newvar s
rome_conditional_constraint(s <= r(1:50)'*x, U2, U3);

rome_maximize(s);   % Note that r is uncertain and depends in z.     
                    % The objective should be interpreted as 
                    % max{y : y<= r(z)'x for z in uncertainty set}

% solve
h.solve;
robust_obj(3) = h.objective;
exp_return(3) = mu(1:50)'*h.eval(x);

rome_end;


load +RomeTest/ConditionalConstraintData.mat

% Error Check:
robust_diff = robust_obj - robust_obj_benchmark;
exp_diff = exp_return - exp_return_benchmark;

% report
disp(sprintf('Case 1, robust err = %0.2f, expected err = %0.2f', robust_diff(1), exp_diff(1)));
disp(sprintf('Case 2, robust err = %0.2f, expected err = %0.2f', robust_diff(2), exp_diff(2)));
disp(sprintf('Case 3, robust err = %0.2f, expected err = %0.2f', robust_diff(3), exp_diff(3)));

% remove all unnecessary variables
clearvars -except ROME_ENV

